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Abstract 

Background: Paired survival data are often used in clinical research to assess the prognostic effect of an exposure. 
Matching generates correlated censored data expecting that the paired subjects just differ from the exposure. 
Creating pairs when the exposure is an event occurring over time could be tricky. We applied a commonly used 
method, Method 1 , which creates pairs a posteriori and propose an alternative method, Method 2, which creates pairs 
in "real-time". We used two semi-parametric models devoted to correlated censored data to estimate the average 
effect of the exposure HR(t): the Holt and Prentice {HP), and the Lee Wei and Amato {LWA) models. Contrary to the 
HP, the LWA allowed adjustment for the matching covariates {LWA a ) and for an interaction {LWAj) between exposure 
and covariates (assimilated to prognostic profiles). The aim of our study was to compare the performances of each 
model according to the two matching methods. 

Methods: Extensive simulations were conducted. We simulated cohort data sets on which we applied the two 
matching methods, the HP and the LWA. We used our conclusions to assess the prognostic effect of subsequent 
pregnancy after treatment for breast cancer in a female cohort treated and followed up in eight french hospitals. 

Results: In terms of bias and RMSE, Method 2 performed better than Method 1 in designing the pairs, and LWA a was 
the best model for all the situations except when there was an interaction between exposure and covariates, for 
which LWAj was more appropriate. On our real data set, we found opposite effects of pregnancy according to the six 
prognostic profiles, but none were statistically significant. We probably lacked statistical power or reached the limits of 
our approach. The pairs' censoring options chosen for combination Method 2 - LWA had to be compared with others. 

Conclusions: Correlated censored data designing by Method 2 seemed to be the most pertinent method to create 
pairs, when the criterion, which characterized the pair, was an exposure occurring over time. In such a setting, the 
LWA was the most appropriate model. 
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Background 

Five percents of breast cancers occur among women 
before the age of 40 and in fewer than 2% of women 
under 35 [1]. Physicians are more and more faced with 
issues of post-treatment pregnancy, such as the optimal 
time before considering conception. In physiopathologi- 
cal terms, a pregnancy after breast cancer is not advisable, 
especially in patients with positive hormonal receptors 
[2,3]. However, many observational studies report that 
pregnancy does not have any impact on the evolution of 
breast cancer [4-6], and could even have a long-term pro- 
tective effect [2,6-10]. This phenomenon can be due to 
the so-called "healthy mother bias" [7], i.e. only women 
who are in good health will undertake a pregnancy after 
treatment for breast cancer. 

Pregnancy is therefore related to the prognostic con- 
dition of the patient, and its occurrence might change 
the prognostic effect of some factors, for instance, the 
hormonal receptors. The prognosis of positive hormonal 
receptors could be qualitatively (interaction) different 
before and after pregnancy {i.e. interaction between hor- 
monal receptor and pregnancy). 

In an attempt to control for this "healthy mother effect" 
various statistical approaches have been used. Some 
authors used the standard Cox model [11] with preg- 
nancy considered as a covariate whose value depends on 
time, and adjusting for the known prognostic factors at 
the diagnosis of cancer, reflecting the gravity of the dis- 
ease [2,6,8,9]. In spite of this adjustment for disease grade 
at diagnosis, pregnancy still remained a significant good 
prognostic factor, indicating a long-term protective effect 
of pregnancy, which is difficult to explain and to interpret. 
Others tackled the problem from the angle of an "illness- 
death" model [12-14] which makes it possible to describe 
the natural history of the disease, taking into account the 
prognostic profiles of the patients who had or did not 
have a pregnancy [15]. This model provided better under- 
standing and interpretation of the effect of pregnancy 
by comparing transition probabilities to relapse between 
women who had or did not have a pregnancy. Moreover, 
taking into account a possible interaction between prog- 
nostic profiles and the pregnancy improved the estimation 
of the pregnancy prognostic effect. Overall prognosis was 
not adversely affected by subsequent pregnancy. How- 
ever, to allow adjustment, such a complex model requires 
enough pregnancies and events of interest. 

In both previous approaches, the adjustment relied only 
on known prognostic factors. However, measured prog- 
nostic factors might not be enough to characterize the 
prognostic status in such a disease. Matching subjects 
designs might be helpful to accomplish that. 

Other researchers have conducted paired studies: preg- 
nant and non-pregnant women were matched on the 
main known prognostic factors (hormonal receptor, 



proliferation level, nodal involvement, use of chemother- 
apy, year of diagnosis), and the non-pregnant had to be 
disease-free for as long as the time from diagnosis to preg- 
nancy of the pregnant women [3-5,7,10]. By this matching 
carried out on known and measured factors, one can 
suppose that the subjects of the same pair also share 
non-observable, not observed or not measured factors, 
in addition to the factors of pairing. Thus, this design 
may improve the control of the "healthy mother effect" 
compared to the two approaches presented above. 

However, to our knowledge, in such a case and con- 
trary to the first methods cited previously, researchers 
[3-5,7,10] did not take into account the fact that preg- 
nancy was an event occurring over time. They matched 
the pregnant woman to a non-pregnant one a posteri- 
ori, i.e. at the end of the follow up study, knowing which 
women were pregnant and which were not over the study 
period. They analyzed the data as if these pairs were a 
priori known and created at diagnosis, i.e. at time t = 
0. Moreover, they always used the stratified Holt and 
Prentice semi-parametric model (HP) [16] to estimate 
the pregnancy prognostic effect, whereas other semi- 
parametric models devoted to censored correlated data 
are available such as frailty models [17,18] and marginal 
models [19-22]. Frailty models model the time distribu- 
tion conditionally to a random effect (frailty covariate), 
specific to each pair, and which is not observed. The struc- 
ture of correlation has to be defined. The latter leave the 
nature of dependence among paired failure times com- 
pletely unspecified. 

Non-parametric [23,24] and parametric [25] approaches 
have been developed, but we focus on the semi- 
parametric approach, more specifically on the commonly 
used marginal semi-parametric model. This marginal 
approach was developed by Wei, Lin and Weissfeld [19] to 
analyze subjects with multiple events, and then Lee, Wei 
and Amato [20] adapted it to clustered subjects. 

In this paper, we use the marginal paired proportional 
hazards model of Lee, Wei and Amato (LWA) [20] and the 
Holt and Prentice stratified model [16] (HP). The main 
difference between them lies in the ability of the LWA [20] 
model to adjust for matching covariates and for the pos- 
sible interaction between the covariate and the exposure, 
contrary to the HP model. Mehrotra et al. [26] proposed 
an efficient alternative to the stratified Cox model analysis 
to estimate the exposure effect, which does not require the 
assumption of a common hazard ratio across strata. How- 
ever, that model is not adapted to our particular context of 
a large number of strata, with very small sample size per 
stratum (in our work, a stratum is a pair), thus it will not 
be studied here. 

With HP and LWA models, we considered two differ- 
ent methods to create our pairs: the a posteriori one 
commonly used and described previously, and a new one 
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designing the pairs "in real-time" by taking into account 
the occurrence of the event over time, i.e. the preg- 
nancy, which characterizes the subjects group within the 
pair. 

The goal was to determine the combination between 
matching methods {a posteriori and in "real-time") and 
models (HP and LWA), which is the most efficient in terms 
of bias and Root Mean Square Error (RMSE) to estimate 
and test the pregnancy prognostic effect. In the follow- 
ing, pregnancy is referred to as the "exposure" occurring 
over time and the pairs are composed of an exposed and a 
non-exposed subject. 

In the next section, the two matching methods as 
well as the models of analysis are presented. In Section 
'Simulations', extensive simulations are used to analyze 
the performance of these Method - Model combina- 
tions in term of bias and RMSE of the estimated effect 
related to the exposure occurring over time. In Section 
'Real data application', our findings are applied to real 
data in order to analyze the effect of a subsequent preg- 
nancy on breast cancer evolution [9]. A brief discussion 
concludes the paper. 

Methods 

The subjects are matched on all the known prognostic 
factors represented by vector Z. In the following, index j 
corresponds to the rank of the exposed subject occurring 
at the yth positions order among the n e exposed subjects 
at the end of the study (j = 1, . . . , n e ), and index i corre- 
sponds to the subject among the n subjects of the study 
(i = 1, . . . , n). Let t{ being defined as the follow-up time of 
each subject /, and fe, as the time exposure for a subject /; 
for a subject who is never exposed, we consider fe; = +oo. 
R m (t) represents the subjects at risk of event and non- 
exposed at time t for the method m (m = 1 or 2). Let the 
set of paired subjects being defined as Pj = {/, /(/)}, with 
the subject i of R m (fe ; .) chosen for the yth exposed 
subject. 

The two matching methods 

The two matching methods applied in the context 
where exposure is an event occurring over time, are 
presented. 

Method 1 matches the subjects as follows: all the sub- 
jects i exposed (fe; < +oo) are matched with a subject 
who was not exposed during the whole period of follow- 
up. This method is called an a posteriori method because 
the groups exposed and non-exposed are defined at the 
end of the study and are considered to be known and 
created at time t = 0. According to Method 1, if a 
subject y undergoes an exposure at time fey, then the 
eligible subjects set R\ (fey) of subjects i eligible to be 
matched to y could be written as follows: R\ (fey) = 
{/ z£ j/ti > feyANDfej = +oo }. This approach has been 



commonly used in the literature [3-5,7,10] with exposure 
not being considered as time-dependent. 

As exposure is a time-dependent event, a second 
approach is proposed which designs the pairs "in real- 
time": according to Method 2, if a subject j undergoes an 
exposure at time fey, then the eligible subjects set R2 (fey) 
of subjects i eligible to be matched to y at time fey, could 
be written as follows: R 2 (fey) = {i 7^ j/U > t£j AND 
fe > ^2 (fey) includes subjects at risk that are not 
yet exposed at fey and will be exposed after, and subjects 
who will never be exposed. A pair composed of an exposed 
subject and a non-exposed one who would never undergo 
the exposure is called a "perfect pair", whereas it would be 
an "imperfect pair". As a result, the non-exposed subject 
of an imperfect pair could be matched after exposure with 
another non-exposed subject. 

Method 2 is similar to the one proposed by Lu et al 
in case-cohort studies [22] where the membership in 
the exposed and unexposed groups is the outcome to 
be explained, whereas in our work it is an explanatory 
variable. 

For both methods, exposed and not exposed subjects 
are matched according to the covariate vector Z, and the 
not exposed subject has to be disease-free for as long as 
the time from the starting point (disease diagnosis) to the 
exposure time. 

For both methods, if several non-exposed subjects can 
be matched with an exposed one, the matched non- 
exposed subject is randomly chosen from the set of 
eligible subjects R m (t); if no non-exposed subjects are 
available, the exposed subject cannot be paired and is thus 
excluded from the analysis. 

Note that even if a subject could belong to two differ- 
ent pairs with Method 2, these two pairs are independent 
while they are never at risk at the same time t. 

The statistical models 

In the following, Xi(t) is the instantaneous hazard function 
of outcome to be estimated for pair Pj. It is noted A.,- (£, Z/) 
to specify that the estimation is made on the pair Py, which 
is composed of the exposed subject j and the non-exposed 
subject i matching on Z/. This notation is the same for all 
the models studied, even those where the adjustment for 
Zj is not available. 

For all the models presented below, E[ it) corresponds 
to the time-dependent exposure status and is defined as 
follows: Ei (t) = 0 if t < fe/, and = 1 if t > fe/. 

The pair of subjects is also defined by a time-dependent 
covariate: Pi (t) = j if i e Pj and t e[t£j; U] , or Pi (t) = 0 
otherwise. 

Holt and Prentice stratified Cox model 

Holt and Prentice [16] adapted the standard Cox model 
[11] to analyze matched paired data. 
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The instantaneous hazard function is written for each 
subject i as 

ki (t, Zd = k m it) exp (y (t) Ei it)) (HP) 

^oi(j) if) is a pair-specific baseline hazard function that 
is assumed to be identical for both subjects of pair Pj, con- 
sidered here as strata; it is considered as a nuisance param- 
eter not to be estimated. The exposure effect exp (y (t)) 
is then estimated, considering the between-pair hetero- 
geneity, by allowing the instantaneous baseline hazard to 
be different within each pair. It is assumed to be iden- 
tical across strata (no interaction between the exposure 
and the pairs) and thus to be implicitly common for the 
whole exposed population: exp (y (t)) is defined as the 
population-weighted average of the stratum-specific haz- 
ard ratios. However, if this assumption is incorrect, i.e. 
in the presence of a true (and often undetected) interac- 
tion, using this model leads possibly to a biased and/or less 
powerful analysis [26]. 

Furthermore, with this model, estimation of the expo- 
sure effect cannot be adjusted for a possible interaction 
between the matching factors and the exposure. This 
stratified approach is sensitive to the unit number per 
strata and to the number of strata: the accuracy of the 
regression coefficients decreases for a small number of 
units per strata and/or many numbers of strata [27]. 

This model is implemented in R software [28] through 
the coxph function by including the term " strataiP lit))" 
with the other explanatory covariates. 

Lee, Wei and Amato Cox model 

The marginal LWA model [20] is an alternative to the 
standard Cox model [11] and is written as follows 

ki it, Zd = k 0 it) exp (y (t) Ei (t)) , (LWA U ) 

if the exposure effect is not adjusted for the matching 
covariates vector Z; 

kt it, Zi) = k 0 it) exp (fi'Zi + y it) E t it)) , (LWA a ) 

if the exposure effect is adjusted for the matching covari- 
ates vector Z; 

ki it, Zi) = k 0 it) exp (fi'Zi + y (t) E t (t) + a! it) Z t E t (t)) , 

(LWAi) 

if the exposure effect is adjusted for the matching covari- 
ates vector Z, and for the interaction between Z and the 
exposure. 

For each of these three LWA models, ko (t) is an unspec- 
ified marginal baseline hazard function considered as 
common for all the pairs, so for the whole population. As 
above, it is considered as a nuisance parameter; exp (y (t)) 
is the average time-varying exposure effect as in the HP 
model, but adjusted (LWA a ) or not (LWA U ) for covariates 
Z and for the possible interaction between covariates and 



exposure (LWAi). Like the standard Cox model [11], the 
LWA assumes that all sample subjects are homogeneous 
(all subjects have the same ko (t)) in spite of the possi- 
ble adjustment for covariates (unique difference between 
LWA U and HP). 

This model is implemented in R software through the 
coxph function, by including the term " cluster (P lit))" with 
the other explanatory covariates. 

For both models, the Proportional Hazard Assumption 
(PHA) was evaluated by Harrels test on scaled Schoenfeld 
residues. This test is implemented in R software through 
the cox.zph function. The possible time-dependent effect 
of the exposure was taken into account by time intervals 
chosen a posteriori, and not by a time-specified function. 

Note that the combination HP and Method 1, taking 
the exposure as a time-dependent covariate or not, gave 
exactly the same estimation of HR (t), whereas LWA did 
not. 

Table 1 presents all the models according to the adjust- 
ment or not for covariates. 

In the following, different simulations and analyses were 
performed with R software version 2.13.0. 

Results 

Simulations 
Objective 

The main objective of the simulation study was to assess 
the ability of the HP and LWA models to estimate the 
true effect of exposure HR it), defined by exp (y (t)), in a 
context of matched paired survival data, where the pairs 
were designed according to the two different methods 
described previously. The aim was to establish the most 
efficient Method - Model combination. 

Dataset 

Simulation of cohort data - Procedures and scenar- 
ios chosen. All the details of the cohort data simulation 
and the procedures and scenarios chosen are given in 
Appendix A. 

We simulated the cohort data referring to an "illness- 
death" model with transition intensities ku it), ku (t) and 
k23 it) (Figure 1). The parameter of interest HR(t) corre- 
sponded to the ratio A.23 (t) /ku (t). The average HR(t) is 
obtained from an exact formula involving the averages of 
A. 13 (t) and A.23 (f) which are computed through a numer- 
ical approximation (transformation of the time from con- 
tinuous to discrete values) (See the Appendix B). The 
average HR(t) adjusted for the different covariates was 
estimated empirically: its estimation was obtained using 
large size samples to guarantee good precision. Moreover, 
note that the larger the ratio ku it) /ku it), the larger the 
number of exposures in the simulated cohort. 

Each subject was characterized by a prognostic profile 
through vector Z, which corresponded to three dummy 
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Table 1 Z, ) estimations with the 7/P and ZWA models are associated to the HR(i) to be estimated according to the 
adjustment for covariates 



Exposure covariate 



Models 



HP 



LWA 



HR(t) 



Time-dependent 

and no time-dependent effect 

With adjustment for Z 

and with interaction 



k m (t) exp (yEj(t)) 



X 0 (t) exp (yEj(t)) 
k 0 (t) exp (y£/(r) + /J'Z/) 
A. 0 (t) exp (/£,- (t) + + «' (t) Z/£, (t)) 



HK(t) 
HRi(t) 



Time-dependent 

and time-dependent effect 

With adjustment for Z 

and with interaction 



A.o/0) (0 ex P (X (0 5(0) 



A. 0 (t) exp (y (r)£/(r)) 
A. 0 (t) exp (y (Of/CO + ZJ'Z/) 
k 0 (t) exp (y (t) E/ (t) + fi'Zj + «' (t) Z/f/ (0) 



H/?(f) 
H/?/(t) 



A./(f, Z,) is the instantaneous risk function of event for the pair Pj, which is composed of the exposed subject j and the non-exposed subject / matching on Z,. 



covariates (k = 1,2 or 3). At time t = 0, there were 
2 3 = 8 possible profiles of Z factors. 

The simulation model included (i) the choice of an 
instantaneous baseline risk function k uv (£, Z) for each 
of the three transitions u -> v (Table 2), (ii) the 
choice of the Z effects, exp f° r eacn transition, i.e. 

k uv (t, Z) = X 0>uv (t) exp (PuviZi + ^ Mv2 Z 2 + PUV3Z3) and 
(iii) the choice for the censoring proportion. 

For (i), an instantaneous average risk function 
k uv (t, Z = Z) for each of the three transitions was sim- 
ulated. Table 2 displays the k uv (t, Z = Z) distributions 
of each transition used for each of the five different 
configurations oiHR (t). 

For (ii), ten different /3 uv k scenarios considered as plau- 
sible fiuvk clinical values [9,15], were performed. Given the 
five configurations chosen for HR(t) and these ten p uv fr 
scenarios, 50 different situations were obtained. 

Finally, for (iii), these previous 50 situations were 
first performed without censoring. Two levels of inde- 
pendent uniform censoring were implemented only to 



the following p uvk scenario: f$ 12 = (-0.2, -0.4, -0.8), 
P 13 = ~p 12 and P23 = Pn> an d they were applied to 
each of the five configurations of HR (t). This yielded to 
10 more situations. 

For each of the 60 situations, 1000 different data sets 
were generated with a sample size of 2000 subjects. At 
t = 0, these 2000 subjects were allocated to eight Z pro- 
files. At t > 0, the 250 subjects of the 8 different profiles 
will be divided up in the three transitions and will change 
over time according to the five HR (t) configurations. 

All theoretical values of HR (t) were calculated on the 
simulated cohort data. They were computed in the over- 
all correlated censored data and inside each sample of the 
Z profile. The average HR (t) was calculated without and 
with adjustment for the matching covariates, i.e. HR (t) 
and HR a (t) respectively. HRi(t) is the one estimated in 
each Z profile. 

From the correlated censored data, HP and LWA U mod- 
els are both assumed to give an average HR(t) i.e. HR (t) 
(considering different assumptions), so they are the only 




State 2 

Pregnancy 



State 1 

Breast cancer 




State 3 




Outcome 


diagnosis 





Figure 1 "Illness-death" model. "Illness-death" model with three transition intensities k uv (t). 
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Table 2 Survival functions applied to simulate each transition and each selected configuration 





Constant H/?(f)* 


Increasing HR(t) 


Decreasing Hfl (f) 


Increasing then 
decreasing HR (f) 


Transition 1 ->> 2 


Exponential^ 


Exponential (1) 


Exponential^ 1 ) 


Exponential^ 1 ) 




(A) 

(0.0050) 


(A) 

(0.0075) 


(A) 

(0.0025) 


(A) 

(0.0020) 


Transition 1 3 


Weibull (2) 


Weibull (2) 


Loglogistic (3 ) 


Weibull (2) 




(W) 

(0.0039,1.1881) 


(W) 
(0.0022,1.1881) 


(5.7146,0.2390) 


(A 0 ,y) 
(0.0018,1.1881) 


Transition 2^3 


Weibull (2) 


Weibull (2) 


Loglogistic (3) 


Loglogistic (3) 




(W) 

(0.0039,1.1881) 


(W) 

(0.0028, 1 .5439) 


(5.6778,0.2463) 


(M,o-) 
(5.9858,0.4971) 



(DS(f) = exp(-At); ^S(f) = exp(-(A 0 t) y ); (3) $(0 =[1 +exp(^^)]" 1 . 

*We simulated the same functions and parameters for the second Constant HR(t) except for Transition 2^3 where k 0 = 0.0069. 



two models which could be compared. Fitting of model 
LWA a makes it possible to estimate an average HR(t) i.e. 
HR a (i), while LWA{ is assumed to give HR[{t) for each Z 
profile (Table 1). 

As the exposure effect was considered to change over 
time for three of the five configurations, its estimation was 
assessed by time interval specified a posteriori. 

Matching methods - Creation of censored correlated 
data from cohort data. For each data set, the two match- 
ing methods presented in Section 'Methods' were applied. 

According to Methods 1 and 2, the two subjects in each 
pair were matched on the three covariates Z/ c , and the 
non-exposed subject had to be disease-free for as long as 
the time from t = 0 to exposure time of the exposed 
subject. 

Then from the 2000 subjects simulated in cohort data 
sets and equally allocated to the 8 Z profiles, a number of 
pairs smaller or equal to the number of subjects in State 
2 (i.e. pregnancy) were obtained. This latter depended on 
the situation simulated, resulting from the HR (t) config- 
uration, the /3 uv k scenario and the censoring percent. 

Statistical criteria used to compare the perfor- 
mances of the different estimators. To estimate a time- 
dependent effect, the time interval [0 — £ max ] was divided 
into L time intervals // defined a priori, according to the 
HR(t) configuration, and written as follows: 



uq = 0 < a\ < . . . < ai = t n 



and 

// = [ai-i;ai[,l > 1. 

If there was no interaction and no time-dependent 
effect, an estimation y s was obtained that corresponded 
to the estimation of y performed in the s th simulation set 
inside the same HR(t) configuration. 



If there was an interaction and no time-dependent 
effect, an estimation for each of the 8 Z prognostic profiles 
was obtained and expressed as y s + a s Z. 

If there was no interaction and a time-dependent effect, 
an estimation for each of the // time intervals was obtained 
and expressed as y s \. 

If there was an interaction and a time-dependent effect, 
8xL estimations were obtained and expressed as y s \ + 

«>• 

To assess the combination Method - Model to estimate 
HR (t), the bias and the RMSE of the estimations pre- 
sented above were calculated. The 8 Z profiles are indexed 
by w = 1 to 8, and the profile w is noted Z^ w \ 

In the event of with an interaction and a time-dependent 
effect, the bias was estimated for each / and each Z^ over 
the 1000 simulations as follows: 

^ 1000 

b izw = J2 + "si z - (ri + <*i z )) 



S=l 



Yl 



+ a , / Z-(y / + « , / z) 



and 



RMSE lz(w) = Jb 2 aiw) + V(y>i + &',Z), 



where, in each // interval, y\ + o^Z is the mean and 
V i^yi + <£/Z^ the empirical variance of the 1000 parame- 
ters estimated y s \ + a 5 /Z: 



V 



(1000 - 1) 

_ [riff (&+ &>)] 

1000 



Eiooo / ^ y 



2~1 
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To compare the different combinations Method - 
Model, the bias was averaged over the profiles or over 
the time intervals (b*z^) or both (b..). 



1 1 

r(w) = J^2b lz (w) 
L 1=1 

w=l 

i 8 



w=l 



The associated RMSEs are given by 



RMSE 9Z(W) = 



RMSEi. = 



RMSE. 



N 1 /=! 



1 

8E^) + ^(M + «i Z ) 



Note that if the exposure effect is not time-dependent, 
then L = 1. If there is no interaction between the exposure 
and the prognostic profiles, then a = 0. 

/?esu/ts 

One particular situation among the 60 simulated is 
described. Figure 2 displays the increasing then decreas- 
ing HR (t) configuration, for each profile and on average, 
without censoring and with 0 12 = (—0.2, —0.4, —0.8), 
P 13 = — 0 12 and 0 2 3 = —/J13. In this particular sit- 
uation which corresponds to a "healthy effect" because 
of the negative values of pu, Figure 2 shows three 
different overall effects of the exposure: a pejorative 
one in the three better prognostic profiles (PP) (Z e 
{(0, 0, 0) , (0, 0, 1) , (0, 1, 0)}), no effect in the intermediate 
PP (Z' e {(0, 1, 1)}) and a protective effect in the last 
fourPP(Z e {(1,0,0), (1,0,1), (1,1,0), (1,1,1)}). With 
023 = ~Pu> we force an interaction between Z and the 
exposure. Note that in this particular configuration cho- 
sen, where 0 2 3 = ~Pi3> HR it) ~ HR a (t) and their values 
are so close that the difference between them is not visible 
in Figure 2. 

Number of pairs. Inside each profile, the maximum 
number of pairs was determined by the number of 
exposed subjects. With Method 1, this number was also 
limited by the number of "perfect" non-exposed subjects, 
but not with Method 2 since the non-exposed subject set 



was composed of "perfect" and "imperfect" non-exposed 
subjects. The difference between the number of pairs 
from the two methods depended on the number of expo- 
sures: the larger the number, the larger the difference. We 
computed the relative difference (RD) in number of pairs 
between the two methods defined as 



RD - 



Number of pairs with Method 2 — Number of pairs with Method 1 
Number of pairs with Method 1 



whose median was equal to +55% (range, +23% 
to +220%). Figure 3 represents the distribution of the 
number of pairs according to the profiles and to the 
matching methods: the median number with Method 2 
was always larger than or equal to that with Method 1. 

Figures 4 shows the number of subjects pertaining to the 
three possible subjects groups at each time t: the exposed 
subject (solid green line), the non-exposed subject who 
never will undergo the exposure (solid red line) and the 
non-exposed subject who will undergo the exposure 
(solid blue line). The dotted vertical green line represents 
the time of first exposure, i.e. the time of occurrence of 
the first pair; the dotted vertical blue line corresponds to 
the time of the last perfect pairs creation and the dot- 
ted vertical red line corresponds to the time of the last 
imperfect pairs creation. With Method 2, the larger the 
ratio X\2 (t) /A. 13 (£), the larger the number of imperfect 
pairs and thus the greater the probability for an exposed 
subject to belong to an imperfect pair. Table 3 provides 
the proportion of imperfect pairs among the whole pairs 
which was estimated over the 1000 simulated data sets of 
our particular situation. It was equal to 81% in the good 
profile Z = (0,0,0), decreasing to 54%, 44% and 20% in 
the Z profiles (1,1,0), (0,0,1) and (1,1,1), respectively. 
Moreover, the larger the ratio A12 (t) /A43 (£), the faster 
the pairs creation stopped. After the last dotted line (blue 
or red, depending on the Z profile), the exposed subjects 
are no longer able to be matched with a non-exposed one, 
because they are no longer available subjects; the pairs 
creation stopped at a time that gradually increased from 
Z = (0,0,0) to Z = (1, 1, 1). For instance, this is illus- 
trated for profile Z = (0, 0, 0) in Figure 4A: at the time 
of first exposure (dotted vertical green line), the exposed 
subject was more likely to be matched with an imperfect 
subject than to a perfect non-exposed one. This set of 
non-exposed subjects decreased over time because each 
of them was matched with an exposed subject until the 
dotted blue line, when no more non-exposed subjects 
were available, while a new exposed subject, belonging 
before to a pair as a non-exposed one, appeared. This 
set of exposed subjects increased and was not able to be 
matched because there were no longer any non-exposed 
subjects. For a same level of A.12 if) A13 (t) ratio, RD was 
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Figure 2 HR (t) configuration chosen. Increasing then decreasing HR (t) configuration, for each Z profile and on average, without censoring and 
with /? 12 = (-0.2, -0.4, -0.8), /? 13 = -/? 12 an ^ $23 = -/? 13 .This figure displays the theoretical estimations of HR(t) called "mean", HR a (t) called 
"adjusted mean" and /-//?,■ (f) in the eight prognostic profiles. The profile Z = (0, 0, 0) at time t = 0 is the profile with the better prognosis; the profile 
Z' = (1,1,1) has the worse prognosis, and the 6 others an intermediate prognosis. In this particular configuration chosen, where /? 2 3 = 
HR (t) ~ HR a (t) and their values are so close that the difference between them is not visible in this figure. 



Numbers of pairs 
according to the profiles and the pairs design's methods 



■ M 2 





(0, 0, 0) (0,0,1) (0,1,0) (0,1,1) (1,0,0) (1,0,1) (1,1,0) (1,1,1) 

Profiles 

Figure 3 Number of pairs. Distribution of the number of pairs according to the profiles and to the matching methods M] and Mj. Results obtained 
with the increasing then decreasing HR(t) configuration, without censoring and with fi u = (-0.2, -0.4, — 0.8), ft ]3 = -/? 12 and /? 23 = — $13. 



Savignoni etal. BMC Medical Research Methodology 201 4, 14:83 
http://www.biomedcentral.eom/1 471-2288/1 4/83 



Page 9 of 18 



CD 
-Q 
E 
13 



O 
CD 

13 
C/) 



Profile (0, 0, 0) 




Perfect non-exposed 

Imperfect non-exposed 

Exposed 



200 400 600 

Time since diagnosis 



800 



1000 



Profile (0, 0, 1) 



j bjects number 

D 100 150 200 

i i i 


D 
D 






- Perfect non-exposed 

- Imperfect non-exposed 

- Exposed 










CO 10 " 










o - 







200 400 600 

Time since diagnosis 



800 



1000 



Profile (1, 1,0) 



Perfect non-exposed 

Imperfect non-exposed 

Exposed 



CD 
.Q 

E 

13 
C 
(f) 
O 
CD 

ZT 
13 

(f) 




CD 
.Q 

E 
13 
C 

cn 
o 

CD 

ZT 

13 

V) 



Profile (1, 1,1) 



Perfect non-exposed 

Imperfect non-exposed 

Exposed 




200 



400 



600 



800 



1000 



Time since diagnosis Time since diagnosis 

Figure 4 Number of subjects in the three possible groups. Number of subjects pertaining to the three possible subjects groups at each time t: 
the exposed subject (solid green line), the non-exposed subject who never will undergo the exposure (solid red line) and the non-exposed subject 
who will undergo the exposure (solid blue line). The dotted vertical green line represents the time of first exposure, i.e. the time of occurrence of the 
first pair; the dotted vertical blue line corresponds the time of the last perfect pair's creation; and the dotted vertical red line corresponds the time of 
the last imperfect pair's creation. Results obtained with the increasing then decreasing HR(f) configuration, without censoring and with 
Pn = (-0.2, -0.4, -0.8), 0 13 = -p v and = -0 13 . 



larger depending on the censoring percent: the higher the 
censoring percent, the smaller the RD. 

Bias and RMSE. Figures 5 displays the distribution of 
the bias and the RMSE of the exposure effect estimator, 
according to our four models. 

With Method 1, the estimation of HR (t) was biased 
with HP and LWA U , but more with HP; the estimation 
of HR a (t) was biased with LWA a . Whatever the model 
applied, HR(t) {HP, LWA U ) and HR a (t) (LWA a ) were 
underestimated, i.e. a poor prognostic exposure effect 
tended to be ignored and a no exposure effect tended to 
become a protective one. In each Z profile, LWAi always 
largely underestimated HRi(t) (Figures 5). The bias b. m is 
equal to -0.77, -0.61, -0.58 and -0.52 under HP, LWA U1 
LWA a and LWAi models, respectively; and the RMSE.. is 
equal to 0.88, 0.65, 0.62 and 1.06 under HP, LWA U , LWA a 
and LWAi models, respectively (Figures 5). The bias and 
RMSE are given by PP in Table 3. 



With Method 2 where HR (t) ~ HR a (£), LWA U and 
LWA a gave very similar values, whereas HP gave a smaller 
estimation of HR (t) than LWA U . All these three models 
are very close to the true value of the effect of exposure, 
but none gave an unbiased estimation: this estimation was 
slightly over- or under-estimated according to the time, 
but much less than with Method 1. HRi(t) was overesti- 
mated in some Z profiles: HRi(t) was overestimated in 
the three PP, Z' e {(0, 0, 0) , (1, 0, 0) , (0, 1, 0)}, where the 
ratio A.12 (t) (t) was large and then the proportion of 
imperfect pairs in the set of pairs to be analyzed, was also 
large (Figures 5). The mean of the time-dependent bias b„ 
was equal to -0.10, 0.14, 0.15 and 0.04 under HP, LWA U , 
LWA a and LWAi models, respectively; and RMSE %% was 
equal to 0.40, 0.27, 0.27 and 1.05 under HP, LWA U , LWA a 
and LWAi models, respectively (Figures 5). The biases 
and RMSE are quite acceptable and much smaller than 
with Method 1. Table 3 displays the bias and RMSE more 
specifically, with the percent of imperfect pairs according 



Table 3 Simulations results: time-dependent bias b mZ ( W ) with RMSE mZ ( W ) for the LWAi model and the mean of the number of pairs inside each prognostic profile 
and according to the matching methods are presented with the percent of imperfect pairs for Method 2 



Prognostic profiles Z (M/) % of exposure Method 1 Method 2 







b 9Z ( W ) 


RMSE 9liw) 


Mean of pairs number 


b 9Z ( W ) 


RMSE 9Z(W) 


Mean of pairs number 


% of imperfect pairs 


(0, 0, 0) 


0.82 


-0.97 


1.05 


40.26 


0.13 


0.32 


146.88 


0.81 


(0,0,1) 


0.48 


-0.46 


0.59 


85.46 


0.04 


0.33 


109.05 


0.44 


(0,1,0) 


0.67 


-0.70 


0.79 


68.82 


0.08 


0.31 


133.47 


0.64 


(0,1,1) 


0.30 


-0.22 


0.46 


70.87 


-0.03 


0.36 


74.07 


0.26 


(1,0, 0) 


0.76 


-0.83 


0.91 


54.48 


0.10 


0.31 


141.70 


0.73 


(1,0,1) 


0.39 


-0.34 


0.51 


82.27 


0.007 


0.34 


92.61 


0.35 


(1,1,0) 


0.58 


-0.57 


0.67 


79.94 


0.05 


0.31 


123.44 


0.54 


(1,1,1) 


0.23 


-0.09 


0.46 


56.15 


-0.06 


0.39 


56.88 


0.20 



The percent of exposure inside each prognostic profile is also given. 
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Figure 5 Methods and Models. Exposure effect's estimation according to Methods 1 and 2 and related to: (A) The HP model, (B) The LWA U model, 
(C) The LWA a model and (D) The LWAj model. The solid thick lines red, blue and green (Figure A, B and C) respectively, represent the theoretical 
average of y(t) and the solid thick orange lines (Figure D) represent the theoretical y(t) + c/(t)Z inside each profile Z. The fine lines red, blue and 
green (Figure A, B and C) respectively, represent y(t), estimated according to Methods 1 (fine dotted lines) and 2 (solid dotted lines); and the fine 
red lines (Figure D) represent the observed y(t) + a (t)Z inside each profile. The average bias and RMSE are given for each model: b 99l RMSEi. 
and RMSE.,. These results are obtained with the increasing then decreasing HR(t) configuration, without censoring and with 
012 = (-0.2, -0.4, -0.8), 0! 3 = -p u and = -0 13 . 



to the PP: the larger the percentage of imperfect pairs, the 
larger the bias and RMSE. 

All the conclusions displayed with Method 1 were the 
same for all five configurations and all the p uv k triplet val- 
ues. The biases oiHR (t), HR a (t) and HR( (t) were always 
huge (data not shown) but more or less followed the con- 
figuration and the f} uv k scenarios. All the conclusions with 
Method 2 were valid for all configurations, and for all /3 uv k 
triplet values. In the configurations without interaction, 
i.e. where P 2 3 = Pi3> HP an d LWA U models were more 
appropriate than LWA a and LWAi to estimate HR (t) in 
terms of bias and RMSE, given that LWA U was much less 
biased than HP in most of the scenarios. In some con- 
figurations where the proportion of the profile with the 
smallest HRi (t) was the most highly represented profile 
leading to a low HR (£), HP was better than LWA U . LWA a 



was the only model for estimating HR a (t) and led to very 
slightly biased estimations of HR a (t) (data not shown). In 
the configurations with interaction, i.e. where fi 23 7^ Pi3> 
the LWAi model was the only appropriate one. However, 
in some of these configurations, LWAi slightly biased the 
extreme profiles and not the intermediate ones (data not 
shown). Over the 10 situations with censoring, the cen- 
soring percent ranged from 9% to 48%. Censoring did not 
change any of the previous conclusions. 

Overall, in terms of bias and RMSE, Method 2 per- 
formed better than Method 1 to design the pairs, and 
LWA a was the best model for all the situations except 
when there was an interaction between the covariates and 
the exposure (fi 23 7^ Pw)> f° r which LWAi was more 
appropriate, even if the estimations with HRi(t) were not 
uniformly unbiased. 
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Real data application 
Data 

Our data retrospectively included 870 women treated for 
breast cancer between January 1990 and December 1999, 
and diagnosed before 35 years of age. Information on 
patients' status was collected at the end of the year 2004 
and the median follow-up was 87 months (range, 7 to 
166) at this date. Tumor, treatment and disease evolu- 
tion analyses are available in the paper of Largillier et al. 
[9] . One of the goals of the data analysis was to compare 
disease-free survival between pregnant and non-pregnant 
patients. The protocol was submitted to the appropriate 
French authorities supervising individual computerized 
data files (Commission Nationale Informatique et Liberte 
[CNIL]), and obtained ethical approval from the Institut 
Curie ethics committee (Comite des Etudes en Recherche 
Clinique [CERC]). The GETNA Working Group, which 
was responsible for the data conception, design and acqui- 
sition, allowed us access to these real data. 

We created pairs of pregnant and non-pregnant women 
using the two matching methods. Both women were 
matched according to their cancer gravity level using: 
Scarff-Bloom-Richardson grade (SBR, related to cell pro- 
liferation level), pathological node involvement and hor- 
monal receptor status (RH). They were also matched on 
treatment: hormonotherapy prescribed or not to RH+ 
patients, and chemotherapy administered or not before 
and/or after surgery. Within each pair, the non-pregnant 
woman had to be disease-free for as long as the time 
from diagnosis to pregnancy of the pregnant woman. We 
sought to estimate the effect of subsequent pregnancy 
occurring over time on breast cancer evolution. 

According to the results obtained on the cohort data 
[15] and to the known breast cancer clinical prognostic 
factors in the literature, patients with low cell proliferation 
level (SBR I or II) and no node involvement were consid- 
ered to have a good prognostic profile and likely to plan to 
be pregnant ("healthy mother effect")- According to bio- 
logical assumptions and need for therapy (hormonal treat- 
ment if RH+), women with negative hormonal receptors 
were also more likely to be pregnant. Thus, even if these 
factors were not significant in the paper by Savignoni et al. 
[15], it seemed relevant to estimate the effect of pregnancy 
according to the RH status (and the treatment associated) 
and according to a clinical Prognostic Profile (cPP); a poor 
cPP was defined as an SBR grade III and/or pathologi- 
cal involved nodes, leading to six prognostic profiles: RH 
negative with a good cPP, RH negative with a poor cPP, 
RH positive not treated with a good cPP, RH positive not 
treated with a poor cPP, RH positive treated with a good 
cPP and RH positive treated with a poor cPP. Then the 
effect of the pregnancy was estimated in the whole pop- 
ulation by adjusting or not for the matching factors and 
with regard to the six prognostic profiles by adjusting for 



an interaction between the pregnancy, and the RH status 
(associated with treatment) and the cPP respectively. The 
effect of pregnancy was not adjusted for chemotherapy 
treatment. 

Results 

In view of our simulations, Method 2 was the matching 
method to apply on the cohort data in order to create 
correlated censored data that would give a more accurate 
estimate of the exposure effect. To be able to compare 
our results with previous findings [3-5,7,10], we used the 
HP model on correlated censored data designed from 
Method 1. 

First, we applied the Method 1 — HP combination as 
proposed in the literature [3-5,7,10]. Secondly, we applied 
the HP and LWA with pregnancy and pair as time- 
dependent covariates on the paired survival data created 
with Method 2. With LWA, HR (t) estimation was carried 
out (i) without adjusting for matching covariates (LWA U 
estimates HR (t)), (ii) by adjusting for all the matching 
covariates but without any interaction (LWA a estimates 
HR a (£)) and (iii) by adjusting for all the matching covari- 
ates and with an interaction between pregnancy and 
matching covariates (LWAi estimates HR[(t)). The latter 
could be applied in the event of real or assumed biological 
and clinical interactions. We tested the PHA with Harrels 
test on the two correlated censored data sets and with 
each model. The PHA was verified and we did not add any 
time-effect for pregnancy in the models. HR(t) was com- 
pared to 1 by the Walds test with the appropriate variance 
according to the models. 

In the cohort data, only 668 patients presented no miss- 
ing data for the covariates of interest. Among them, 68 
experienced a subsequent pregnancy. We obtained the 
maximum number of pairs available with Methods 1 and 
2, i.e. 68 pairs: all pregnant women were then matched. 
The 68 pairs were not exactly the same between Meth- 
ods 1 and 2. In Method 2 five pairs were imperfect, which 
represented a low proportion (7.4%). Among the 68 pairs 
obtained with Method 1, 32 women experienced an event 
(progression or death): 16 in the group of patients who 
became pregnant after the breast cancer diagnosis and 
16 in the group of patients who did not. Among the 68 
pairs obtained with Method 2, 29 women experienced an 
event: the same 16 in the group of patients who became 
pregnant after the breast cancer diagnosis and 13 in the 
group of patients who did not. No events occurred in the 
imperfect pairs. Only sixteen pairs (23.5%) were common 
between the two matching methods. The number of pairs 
and the number of final events in the pregnancy and non- 
pregnancy groups are given Table 4, according to the 
six profiles and to the matching methods. The number 
of pairs was divided into imperfect and perfect pairs for 
Method 2. 
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Table 4 Real data: number of pairs and final events according to the matching methods and to the six prognostic profiles 



Prognostic profiles 




Method 1 




Method 2 




Number of pairs 


Number of final events 


Number of pairs 


Number of final events 






Pregnancy/nonpregnancy group 


Imperfect/perfect 


Pregnancy/nonpregnancy group 


RH - and goodcPP 


7 


2/2 


1/6 


2/0 


RH - and bddcPP 


16 


4/2 


1/15 


4/2 


RH + not treated and goodcPP 


12 


4/2 


1/11 


4/1 


RH + not treated and badcPP 


19 


5/5 


0/19 


5/7 


RH + treated and goodcPP 


5 


0/1 


2/3 


0/1 


RH + treated and badcPP 


9 


1/4 


0/9 


1/2 


The number of imperfect pairs is also given for Method 2. 



The HP model applied with Method 1 yielded the fol- 
lowing estimations: exp (y) = 0.90, Cl^% [0.36 — 2.21], 
not significant. Kranick et al. [5] and Veletgas et al. [4] 
concluded likewise as did Azim et al. [3], but the latter 
stratified their analysis on estrogen status. Other authors 
[7,10] found a protective role of pregnancy. All these 
authors used the HP model with Method 1 adjusted for the 
other known prognostic factors. In view of our simulation 
results, we believe that this estimation was underesti- 
mated and that the models should be used on censored 
correlated data created with Method 2. Regarding the lat- 
ter, the HP model yielded a larger value than with Method 
1 but was still not significant (exp (y) = 2.00, Cl^% 
[0.75 — 5.33]). In view of our simulations, we conclude 
that HR(t) was widely underestimated with HP using 
Method 1 and only slightly underestimated with HP using 
Method 2. The increasing value of HR (t) estimated by HP 
between the two matching methods was not surprising. 
Method 2 and the LWA model, with or without adjusting 
for matching covariates (LWA U and LWA a ), yielded simi- 
lar and not statistically significant results: exp (y) = 1.22, 
C/ 95 % [0.61-2.42] and exp (y) = 1.24, C7 95 % [0.62-2.45], 
respectively. The difference in HR (t) values estimated by 
the HP and LWA U models was not statistically significant. 
In such a context, we would like to estimate the proper 
effect of the subsequent pregnancy. It is more relevant 
to estimate HR a (t) than HR (t). Even if LWA t showed no 
significant interaction between the matching covariates 
and the pregnancy, because of the biological and clin- 
ical assumptions, we estimated the effect of pregnancy 
according to the six PP defined above using the LWAi 
model. Table 5 presents exp (y) estimations according to 
the models and their 95% confidence intervals. 

The apparent protective role of pregnancy for the the- 
oretically less favorable prognostic profile "RH+ (treated 
or not) and poor cPP" and its increasing role in risk for the 
theoretically best prognostic profile "RH— and good cPP" 
were surprising. There was no statistical significance, but 
it could be because of a lack of power of the combination 
Method 2 - LWAi. 



Discussion 

In our context of exposure occurring over time, we 
focused on two matching methods: Method 1, commonly 
used in the literature [3-5,7,10], and Method 2, our new 
approach. Method 1 composes the pairs in an a poste- 
riori way where the pairs are in fact considered to be 
known at diagnostic time t = 0. Method 2 creates pairs 
in "real-time". Such matching designs create independence 
between pairs but dependence between the subjects of 
the same pair, and specific analytical methods exist for 
such a situation of correlated censored data. In our work 
we studied two semi-parametric models allowing for the 
stratified Holt and Prentice model [16] (HP) and the Lee, 
Wei and Amato model [20] (LWA). To estimate the aver- 
age exposure effect HR (t) and unlike the LWA, the HP 
did not make it possible to take into account either the 
matching covariates or the possible interaction between 
the matching covariates and the exposure. 

The aim of this study was to analyze the HP and LWA 
models using the two matching methods in order to pro- 
pose the most efficient Method - Model combination to 
estimate and test the prognostic exposure effect HR (t) 
estimated through the models by exp (y (t)). 

In view of our simulations, the relative difference in 
the number of pairs between Method 1 and Method 2 
depends on the ratio (t) /A43 (t) and on the censoring 
percent. Compared to Method 1, the number of pairs was 
equal or larger with Method 2. In terms of bias and RMSE, 
Method 2 is more relevant than Method 1 to design the 
pairs, and LWA a is the best model for all the situations 
except when there is an interaction between the covari- 
ates and the exposure (fi 2 3 7^ Pw)> f° r which LWAi is 
more appropriate even if the HRi(t) estimations are not 
uniformly unbiased. 

In our sample data, we applied Method 1 and HP to 
compare our results with those in the literature [3-5,7,10]. 
According to the simulation results, we applied Method 
2 with HP and LWA. With both matching methods, we 
obtained an equal number of pairs (the maximum avail- 
able) but not the same ones. HP used with Method 1 
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Table 5 exp(y) estimations and their 95% confidence 
interval according to the HP and LWA models and to the 
matching methods, with the p-value of the Wald's Test 

Models Method 1 Method 2 

exp (y) = 0.90 exp (y) = 2.00 

y I j q , . C/ 95 = [0.36 - 2.21] C/ 95 = [0.75 - 5.33] 
Holt and Prentice model 

p - value = 0.82 p - value = 0.1 7 

se = 0.46 se = 0.50 



The estimation of the pregnancy effect standard error (se) by the HP and LWA 
models are given. 



gave the smallest estimation of HR (t). HR (t) estimations 
were not statistically different from 1 with HP, LWA U and 
LWA a . According to biological and clinical hypothesis, we 
assumed that the effect of exposure was different accord- 
ing to the six prognostic profiles so we used the LWAi 
model However, we did not conclude that pregnancy had 
neither a protective effect nor an adverse effect on pro- 
gression disease, even for women with positive hormonal 
receptors and poor clinical profile at diagnosis. Estimat- 
ing HR (t) by prognostic profiles using the combination 
Method 2 - LWAi, showed a different and opposite effect 
of exposure in the six health profiles, but none was sta- 
tistically significant. The sample size of the profiles was 
small and a few events occurred inside each profile, prob- 
ably leading to a lack of power of the combination Method 
2 - LWAi. 

Adjusting for the matching covariates and for their pos- 
sible interaction with exposure could be very interesting 
to interpret the effect of exposure HR (t). However, it 
could be more difficult with real data, especially in the 
event of multiple prognostic profiles. The sample had to 
be large enough to enable us to assess the performances 
of the models according to the two matching methods, 
especially LWA a and LWAi, which required a large num- 
ber of pairs. The erroneous estimation with LWAi in some 
extreme profiles could partially be explained as follows: 
in the "imperfect" pairs, we artificially prevent the non- 
exposed subject from having a final event, because she will 
first become an exposed subject of another pair. Maybe 
she will undergo the final event but as an exposed sub- 
ject and in another pair. Then, when the proportion of 
imperfect pairs is large, HRi(t) could be overestimated. 
Moreover the larger the ratio A. 12 (t) A. 13 (t), the larger the 
number of imperfect pairs, and the higher the overestima- 
tion oiHRiit) and the faster its appearance. Owing to the 
values of the triplets chosen in the particular config- 
uration presented above, the subjects of the three profiles 
Z £ {(0, 0, 0) , (1, 0, 0) , (0, 1, 0)} experienced more expo- 
sure before any events, and simultaneously fewer events 
(before or after exposure), than the other PP. In these 
three profiles, the number of events in the exposed and 
in the non-exposed individuals is quite low, suggesting a 
possible lack of accuracy in the estimation of HR (t). All 
the models needed enough pairs and enough final events 
to fit: the larger the number of imperfect pairs and the 
smaller the number of events, the worse the accuracy of 
the estimation oiHRiit). In a future study, we will explore 
more deeply this problem of bias related to the percent of 
imperfect pairs in the sample. 

Technically, whatever the matching method used, HP 
censored the pair when one of its subjects was censored or 
experienced the final event, whereas the LWA did not. The 
latter leaves each subject to be followed up to her censor- 
ing or final event, whatever the outcome of the matched 



Lee, Wei and Amato model 



LWA U 



exp {y) = 1.22 
C/ 95 = [0.61 - 2.42] 
p — value — 0.58 
se = 0.35 



exp()/) = 1.24 

LWA a C/95 = [0.62 - 2.45] 

p - value = 0.54 
se = 0.35 

exp (y) = 5.70 

LWAi ■ RH- and good cPP C/ 95 = [0.80 - 40.50] 

p - value = 0.08 
se = 0.98 

exp (y) = 1.77 

LWAj : RH- and poor cPP C/95 = [0.34 - 9.24] 

p - value = 0.49 
se = 0.83 

exp (y) = 2.87 

LWAj : RH+ not treated C/ 95 = [0.63 - 1 3.1 0] 

and good cPP 

p - value = 0.17 
se = 0.76 

exp (y) = 0.89 

LWAj : RH+ not treated C/95 = [0.30 - 2.69] 

and poor cPP , „ on 

K p - value = 0.83 

se = 0.55 
exp (y) = 0.95 

LWAj : RH+ treated and C/95 = [0.07 - 1 3.66] 

g00dcPP p-value = 0.97 

se = 1 .33 

exp (y) = 0.30 

LWAj : RH+ treated and C/95 = [0.04 - 2.42] 

poorcPP , 

p - value = 0.25 

se = 1 .05 
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subject. Concerning the management of the imperfect 
pair obtained with Method 2, the non-exposed subject i 
of pair Pj was censored when she became the exposed 
subject i of another pair Pi at the time £g.; by construc- 
tion, as said before, her exposure occurred before any 
other events. For LWA, an alternative would be to cen- 
sor the pair Pj. Another alternative, whatever the model, 
would be to propose another non-exposed (perfect or 
imperfect) subject to the exposed subject of pair Pj, who 
was now single in her pair. This issue deserves more 
investigation. 

Using time intervals to estimate the effect of exposure 
as it changes over time was maybe not the most pertinent 
method, because it required many parameters. Therefore, 
to improve fitting, we intend to apply splines to fit the 
effect of exposure as it changes over time. Indeed, the 
present study sought to compare results of two known 
models devoted to censored correlated data, and the well- 
known frailty model was set aside because it requires the 
structure of correlation within the pairs to be specified. 
The next step would be to compare the HP, LWA and the 
frailty model using Method 2. 

Conclusions 

In conclusion, correlated censored data designing by 
Method 2 seems to be the more pertinent method to cre- 
ate pairs when the criterion which characterizes the pair 
is an exposure occurring over time. It would be inter- 
esting to estimate the proper effect of the subsequent 
pregnancy. It is then more pertinent to estimate HR a (£) 
than HR (£). Thus, LWA a seems to be the best model for all 
the situations, except when there is an interaction between 
the covariates and the exposure, for which LWA{ is more 
appropriate, even if the estimations of ///?;(£) are not uni- 
formly unbiased. LWA a and LWAi gave a more accurate 
and relevant estimation of the effect of exposure in par- 
ticular context, where we can reasonably suppose that the 
latter depends on prognostic profiles. 

Appendix A: Simulation of cohort data - 
Procedures and scenarios chosen 

For each subject four independent times were generated: 
tn> t23> t\3 and C (censoring time), according to the inten- 
sities displayed in Table 2. From these times, 2 times of 
interest for each subject were derived: a time to expo- 
sure tn = t£ and a time to final event t{. Two indicator 
variables were derived: E = 1 if an exposure occurs, 0 
otherwise and A = 1 if a final event occurs, 0 otherwise. 
Four possible quadruplets (££,* t{\ E; A) could be defined as 
follows: 

(C;C;0;0)ifC < min(*i 2 , t 13 ), 
(t 13 ; t 13 ; 0; 1) if t 13 < mm(t 12t C), 
(£12; C; 1; 0) if tn < min(£i 3 , C) and tn + £23 > C, 
ihr, tn + t 23 ; 1; 1) if £12 < £13 and tn + t 23 < C. 



Such a design refers to an "illness-death" model with 
transition intensities A. 12 (t), A43 (£) and A23 (£) (Figure 1). 
All subjects were assumed to be in the initial state (state 
1 or cancer diagnosis in our context) at time £ = 0. 
They could move to the final state (state 3 or disease pro- 
gression) with a transition intensity A. 13 (£). They might 
undergo the intermediate event or exposure (state 2 or 
pregnancy) with intensity A 12 (£), before developing any 
progression with intensity A23 (£). Date of entry into state 
1 was chosen as time of origin for all transitions. Thus 
the parameter of interest HR(t) corresponded to the ratio 
^23 if) A 13 it). However, to compute A23 (£), we took into 
account the left truncation phenomenon: before being at 
risk of an event in the transition 2 -> 3, a subject has 
to wait until its exposure occurs. This delayed entry leads 
the set of subjects at risk in transition 2 -> 3 to increase 
when an exposure occurs and to decrease when an event 
occurs. Thus the average HR(f) is obtained from an exact 
formula involving the averages of A 13 (£) and A23 (£) which 
are computed through a numerical approximation (trans- 
formation of the time from continuous to discrete val- 
ues) (See the Appendix B). The average HR(t) adjusted 
for the different covariates was estimated empirically by 
using large size samples to guarantee good precision. 
Moreover, note that the larger the ratio A 12 (£) A 13 (0> 
the larger the number of exposures in the simulated 
cohort. 

The simulation model included (i) the choice of an 
instantaneous baseline risk function X uv (£, Z) for each of 
the three transitions u —> v, (ii) the choice of the Z effects, 
exp iP uv k)> for each transition and (iii) the choice for the 
censoring proportion. 

For (i), an instantaneous average risk function 
X uv (t, Z = Z) for each of the three transitions was simu- 
lated: either a constant risk using an exponential density 
function^, a monotone risk using a Weibull density 
function^ or an increasing then decreasing risk using a 
loglogistic density function^. Five X uv (t,Z = Z) triplets 
were simulated in order to construct five realistic con- 
figurations of HR (£): two constant, one increasing, one 
decreasing and one increasing then decreasing, where 
HR(t) range values were clinically pertinent (between 
0.5 and 4 in the whole population). Table 2 displays the 
X uv (t, Z = Z) distributions of each transition used for 
each of the five different configurations of HR (£). 

For (ii), different fi uvk values for each of these five 
\ uv (t, Z = Z) triplets were chosen. Negative fin values 
were proposed and set at fi l2 = (—0.2, —0.4, —0.8). 
Only fi 13 and fi 2 3 na d other possible values which 
were the following: (-0.2, -0.4, -0.8), (+0.2, +0.4, +0.8), 
(-0.1,-0.2,-0.4) and (+0.1, +0.2, +0.4) . Ten p uvk sce- 
narios were performed. Given the five configurations cho- 
sen for HR(t) and the ten $ uv k scenarios, 50 different 
situations were obtained. 
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Finally, for (iii), these previous 50 situations were first 
performed without censoring. To minimize simulations 
time, two levels of independent uniform censoring were 
implemented only with the following $ uv k scenario: fi 12 = 
(-0.2, -0.4, -0.8), fi 13 = -fi 12 and 0 23 = fi 12 ; and they 
were applied to each of the five configurations of HR (f). 
This yielded to 10 more situations (five HR (t) configu- 
rations with2 levels of censoring) for that $ uv fc scenario. 
The maximal event time £ max was set at 1000. The first 
uniform distribution for censoring time C was over the 
interval time [0; £ max L an d the second one over [0; 2£ max ]; 
then the maximal censoring time was C max = +oo, £ max 
or 2£ max . The overall censoring level was higher in the first 
censoring distribution but it also depended on the HR (t) 
configuration. In total we had50 situations without cen- 
soring and 10 with censoring (the same five configurations 
with the two levels of censoring). 

For each of the 60 situations, 1000 different data sets 
were generated with a sample size of 2000 subjects. At t = 
0, these 2000 subjects were allocated to the eight Z pro- 
files. At t > 0, the 250 subjects of the 8 different profiles 
are divided up among the three transitions and change 
over time according to the five HR (t) configurations. 

Appendix B: Calcul of the average exposure effect 
not adjusted for the covariates HR(t) 

Density function f 2 3 (t) , Repartition function F23 (t) and 
Survival function S 2 3 if) of transition 2^3 taking into 
account the left truncation 

/23 (t) = f f 23 (t - u)f 12 (u) S 13 (u)du, 
Jo 

rt rt rv 

F23 it) = / f 2 3 (v) dv = / / f 2 3 (v - u)fn (u) Si 3 (u)du, 
Jo Jo Jo 

where u is the exposure occurrence time and t the final 
event time. 

To obtain an useful formula we divide the time in K tiny 



intervals with to < t\ < 



S23 (f*) = fl 1 



< t k < • • • < t K 



f23 (tj) (tj ~ t hl ) 



f = \ \ /o ; /l2 (X) S13 (X)dx - Yli=i (tl ~ t/-l)/23 fe> 



^ ' X fofn (x) S13 (x)dx - Y! l= \ {ti - i/-i)/ 23 m 



k 

n 

7=1 



1 fofn (*) Si 3 (x)dx - (tl ~ ti-i)f23 (ti) 
1=1 

t . j-i 
fofn (x) Su(x)dx -J2(ti~ */-i)/23 (ti) 
1=1 / 



and 



The average instantaneous hazard X (t) is equal to 

_ T,z*>zV)k(t\Z) 

k yt) — ^ — 

CO z (t) = 7T (z)S(t\Z) 

(t) = ^ — 

22z M z (t) 

i(t) = J2^z(t)x(t\z) = j2^z(t)x z (t) 



0) z (t) = 



N z (t) 

Y,N z (t) 



*23(f) 



-10^(^23 (t k )) 
{tk ~ tk-l) 



where X(t\Z) = X z (t) is the average instantaneous haz- 
ard in the profile z at time t, N z (t) is the number of subject 
at risk in the profile z at time t. Then we can write 

" x »<» " • <«»' 

Z 

The computation of N z (t) is the following. 
In transition 1 — ► 3, a subject is at risk if it does not have 
any event and if it is not censored 

Nu z (t)=N 13z S 13z (t)G(t)S 12z (0, 

where 

• Sn z (t) is the survival rate in the transition 1 2 for 
the profile z at time ^, 

• Si3 z (t) is the survival rate in the transitionl — > 3 for 
the profil 2 at time ^, 

• N\3 Z is the number of subject at risk in the profile z at 
time t = 0, 

• G (£) = 1 — G (£) with G(t) is the repartition function 
of the variable C characterizing the censoring time. It 
does not depend on the profile. 

In transition 2^3, taking into account the left trun- 
cation, a subject is at risk of event in this transition, if it 
entered in the transition and if it had neither an event nor 
being censored. 

N 2 3 z could be expressed as 

N 23z (t) = [E 23z (t) - D 23z (t) - C 23z (t)] N 23z 
where 

• E 2 3 z (t): the number of subjects which entered in the 
transition 2 — > 3 in the time interval [0; t] , 

• E>23 z (t): the subjects which undergo the final event in 
the time interval [0; t] , 

• C 2 3 z (t): the censored subject in the time interval 
[0;t], 

• N23 z is the number of subject at risk in the profile z at 
time t = 0. 
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Then 

N 23z it) = [E 23z it) - D 23z (t) - C 23z (t)] N 23z 

~f*G (x) X Uz (x) S Uz (x) S Uz (x) dx 
-foG(x)f*k 12z (v)S Uz (V)Si3 x (v)/ 23z (x-v) dvdx 
-fog( x )fo k i2 z (v)S Uz (v)S Uz (v)S 23z (x - v) dvdx_ 

where 



N 23z 



> feand t e Ii 
otherwise. 



• ^n z (f) is the instantaneous hasard function of 
transition 1 -> 2 for the profil z at time £, 

• g (t) is the density function of the variable C 
characterizing the censored time. It does not depend 
on the profile. 

Appendix C: Empirical estimation of the average 
exposure effect adjusted for the covariates HR a (t) 

Lets consider L intervals // (/ = 1 to L) inside the study 
interval [0; T max ], where all the i subjects are censored at 
time T max : I t = ]ai-i; a{\ : / = 1, • • • , L with 

ao = 0 < a\ < a 2 < • • • < ai = T max . 

Inside each //, the exposure variable £/(£) is as follows: 

^ s , f 1 if ^ 

Then we define: E (t) = X!/=i El if) • 
Our model of simulation obtained from an "illness- 
death" model is as follows: 

X (t, Z,E) = X 0 it) exp (j2 YlEl (f) + ftZi + fcZ 2 
\/=i 

+ p 3 Z 3 + (aiZi + a 2 Z 2 + a 3 Z 3 ) E (t)^j , 

where X(t, Z,E) is the instantaneous risk function of state 
3 occurrence according to the exposure E and the covari- 
ates (/< = 1 to 3). 

From this model, we are able to calculate HRi(t), i.e. the 
HR(t) for each profile (Zi, Z 2) Z 3 ): 

HR i(t) = X(t > Z ' E=1) . 

X(t,Z,E = 0) 

We approximate HR a (t) by 
where the // are obtained from 

X (t, Z,E) = X 0 (t) exp yi £ l W + ^i Z i + + ftZ 3 ^ . 
We approximate by 
M(0 = exp^y / £ / (^, 



where the yi are obtained from 

X it, Z,E) = X 0 it) exp yi £ l V)j • 

For HR a it), the estimation of y/ is adjusted for the /3/ 0 
whereas for HRit) it is not. 

The approximation of HRiit), HR a it) and is espe- 
cially relevant when L and the total number of patients are 
large, in order to have enough exposure and events inside 
the transitions 2^3 and 1—^3. 

The approximations ofHRi it) and///? (t) could be com- 
pared to the theoretical values (See Appendix B), not of 
HR a it). 

To be sure of the approximation of HR a it), we simu- 
lated data from a Cox proportional hazards model with 
a time-dependent covariate representing exposure. As we 
obtained directly the effect of the exposure we were able 
to certify that our approximations method was correct. 
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